Universal scaling of current fluctuations in disordered graphene 
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We analyze the full transport statistics of graphene with smooth disorder at low dopings. First we 
consider the case of one-dimensional (ID) disorder for which the transmission probability distribution 
is given analytically in terms of the graphene-specific mean free path. All current cumulants are 
shown to scale with system parameters (doping, size, disorder strength and correlation length) in 
an identical fashion for large enough systems. In the case of 2D disorder, numerical evidence is 
given for the same kind of identical scaling of all current cumulants, so that the ratio of any two 
such cumulants is universal. Specific universal values are given for the Fano factor, which is smaller 
than the pseudodiffusive value of ballistic graphene {F — 1/3) both for ID {F m 0.243) and 2D 
{F « 0.295) disorders. On the other hand, conductivity in wide samples is shown to grow without 
saturation as \/i and log L with system length L in the ID and 2D cases respectively. 

PACS numbers: 73.20. Jc,73.23.-b 



I. INTRODUCTION 

Stimulated by the striking results of the first ex- 
periments on graphene flakeSfii^ii^ interest in the prob- 
lem of transport in such system has seen extraordi- 
nary growth The particular subject of disorder in 
graphene has been recently the center of numerous stud- 
ies, since many of the known concepts and results for 
transport in disordered normal metals break down for 
the peculiar dispersion relation of undoped graphene. For 
example, simple semiclassical techniques fail;^ quantum 
corrections can show up with opposite sign to the con- 
ventional case (weak antilocalization) unless valley sym- 
metry is broken^iS; such corrections are conspicuously 
absent in the cleanest experimental sampleSfiS the pro- 
posed explanation being an effective time-reversal sym- 
metry breaking due to curvature disorder that preserves 
the valley symmetry in graphene Indeed, the specific 
symmetry properties of disorder turn out to be a cru- 
cial issue for transport in this systemi^iii^ While atomic- 
sized defects are widely thought to be of little importance 
for transport at standard temperatures ji^ii^ii^ii^ smooth 
potentials can have a very visible impact on the conduc- 
tivity (7 of graphene. Such electrostatic potentials can 
arise either from the ubiquitous geometrical corrugation 
observedi^ in most graphene samples;^ or from ineffec- 
tive screening of charges in the environment ii2,ilS, It has 
sometimes been dubbed "charge puddle disorder" due 
to the fact that the gapless nature of graphene makes 
the material respond to such a potential by forming lo- 
cal particle and hole charge accumulations . Recent 
measurements^ estimate the typical size of charge pud- 
dles in the ^ ~ 10 — 30nm range, with a typical potential 
height of ay ^ lOmeV— lOOmeVi ^^i^^ This gives a typical 
dimensionless disorder strength av£,/fivp « 0.2 — 2. 

The effect of smooth disorder on transport in graphene 




[a] ID disorder 



[b] 2D disorder 



FIG. 1: (Color online) Two different disorder realizations con- 
sidered in this work for a graphene sheet (between metallic 
contacts) of length L and width W . (a) ID disorder (no mode 
coupling), (b) model for 2D charge puddle disorder. 



has been recently analyzed theoretically by a number of 
authorS )-^^'^'^i^'^'^^i^^'^''i^^ finding once more striking dif- 
ferences with respect to the well known theory of disor- 
dered metals.— Due to the absence of inter- valley scatter- 
ing and to chirality conservation,^" this kind of disorder 
has the peculiarity of enhancing the conductivity with re- 
spect to the ballistic case, which is at odds with classical 
intuition. A natural question that arises is whether this 
enhancement has an upper bound and whether, as a con- 
sequence, the conductivity of a large graphene sample ex- 
hibits a universal value, as initially claimed,^ - or whether 
it depends on size and system properties, in other words, 
whether the scaling /3(cr) function in smoothly disordered 
graphene exists and has fixed points or not. Some of the 
predictions made so far are conflicting on this point. Us- 
ing a supersymmetric nonlinear sigma model,^ Ref. [2^ 
suggests the existence of a universal minimal conductiv- 
ity, whereas recent numerical simulations find that a in- 
creases in a logarithmic fashion with the system size, Si 
while the beta function exists (single parameter scaling) 
and always remains positivci ^^'^^ One must bear in mind. 
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however, that the conclusions of Ref. |23 rely on a dif- 
fusive limit, so that direct comparison to the numerical 
calculations might not be straightforward. On the other 
hand, a recent experiment"^^ has shed doubts about the 
existence of a universal minimal conductivity in real sam- 
ples. 

In this work, we analyze the effect of smooth puddle 
disorder on the complete transport statistics of graphene 
at low dopings and, in particular, the scaling properties 
of current fluctuations with system parameters such as 
length or disorder strength. First, we study the case 
of one-dimensional (ID) disorder. We numerically com- 
pute transport properties within the transfer matrix for- 
malism. We then derive and solve the single channel 
Dorokhov-Mello-Pereyra-Kumar (DMPK) equation^^ for 
graphene, valid for long samples. The DMPK equation 
implies that a single parameter scaling for conductivity 
holds and /3(cr) = d\oga/d\ogL exists. Both the ana- 
lytical solution and the numerical results agree with zero 
fitting parameters and high accuracy. The most impor- 
tant difference between this result and that of the ID 
disordered metal is that in graphene, we find a channel- 
dependent mean free path that scales quadratically with 
transverse momentum q for small dopings. The conse- 
quence of this in transport for wider-than-long sheets is 
the peculiar scaling properties of the resulting transmis- 
sion probability distribution, which we compute analyti- 
cally. As a consequence of these scaling properties, con- 
ductance and all higher current cumulants scale in the 
same way with system parameters ( "universal scaling" ) , 
such as system length L or disorder strength. In partic- 
ular, conductivity is found to grow without saturation as 
VX with system length, as opposed to the constant con- 
ductivity of ballistic graphene. No sign of localization is 
obtained as expected. The ID Fano factor is found to sat- 
urate, both numerically and analytically, to F = 0.243. 

In the case of two-dimensional (2D) disorder analytical 
progress is difficult. We numerically calculate the trans- 
mission probability distribution, finding that the type of 
scaling features of the ID case are also present in two 
dimenions, albeit with a modified scaling law. Conduc- 
tivity in wider-than-long^D samples is found to scale as 
a oc logi, like in Ref. 24, without localization, in con- 
trast to the initial suggestions.^^ All higher cumulants 
once more exhibit an identical scaling law as the con- 
ductance, leading to truly universal ratios of any pair of 
cumulants, e.g., F = 0.295. 

The layout of this work is as follows. We first give an 
overview of the transfer matrix method employed in this 
work in Sec. |TT1 Then, we study the case of ID disorder 
(Sec. IIII[) . and give a full analytical solution to its trans- 
port statistics at low dopings. The same scaling is found 
explicitly for all current cumulants in this case. Tech- 
nical details about the ID calculation can be found in 
Appendix [B] while a derivation of the DMPK equation 
in graphene can be found in Appendix [Xl Motivated by 
the ID results, we numerically explore the case of proper 
2D disorder in Sec. IIVI and describe the evidence for uni- 



versal scaling of current cumulants also in two dimenions. 
Finally, we conclude in Sec. |Vl 



II. TRANSFER MATRIX METHOD 

Graphene with smooth disorder, which therefore does 
not couple valleys, can be modeled by a single flavor 2D 
Dirac Hamiltonian 



H = vpc ■ p + V{x, y) — e. 



(1) 



where Pauli matrices cr act on the pseudospin space, e 
is the doping, vp ~ lO^m/s is the carrier velocity, p is 
the momentum operator with respect to the Dirac point 
and V{x,y) is a disorder realization (see Fig. [T]). The 
two valleys and two real spins amount to four degenerate 
transport channels in this approximation. In the follow- 
ing, all energy scales such as doping e or potential V will 
be given in units such that hvp = 1 for compactness, so 
that e will in fact stand for e/{hvF), and so on. 

In Refs. [2^ and Is^, Titov and Cheianov and Fal'ko 
derived a differential equation for the transfer matrix 
Tqqi (L) describing the propagation of such Dirac fermions 
with transverse momenta q = 2Trn/W and energy e 
through a graphene sheet of width W and length L un- 
der a given realization of disorder V{x,y). In the case 
of large Fermi wavelength mismatch at the contacts, it 
takes the form^^ 



dx 

Mq'q{x) 
Vq'q{x) 



= y^^Mq'k{x)Tkq{x), 



(2) 



CTx ^qSq'q + ia^ <S) {iSq'q - Vq'q{x)) , (3) 



1 

w 



d2/F(a;,y)e-*(«'-«)^ 



(4) 



The initial condition is T(0) = 1, and dq^qi denotes a 
Kronecker delta. The implicit assumptions in the above 
equations are that the disorder does not couple valleys 
(i.e., it is sufficiently smooth), that metallic contacts 
that connect graphene to the reservoirs are appropri- 
ately modeled by infinitely doped graphene (e — > cx) for 
|a;| > L/2), and that the boundary conditions for the 
transmission modes can be chosen periodic in the trans- 
verse direction. The latter assumption is rigorously valid 
for W ^ L, for which possible boundary- induced pseu- 
dospin precession and valley mixing effects can be ig- 
nored. 

The connection between the transfer and the scattering 
matrix is given by 



T{L) = 



t' 



f't 



where t (r) are transmission (reflection) amplitude ma- 
trices. This allows one to compute eigenvalues T„ of the 
transmission matrix tH, which are conveniently recast 
into parameters A„ through 



sech^Arj 



(5) 
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FIG. 2: (Color online) Normalized transmission distribution 
p(A)/AiD for various ID disordered undoped graphene sys- 
tems. Good agreement is found with Eq. (|13p . plotted in 
red. Inset: comparison of the distribution [Eq. (|8}] to the 
numerical distribution P{Xq) for two different values of the 
transverse momentum q. No fitting parameters were used in 
either case. 



In the context of transport through random disorder 
V{x,y), one is interested in the probability distribution 
of the transfer matrix or, more precisely, the distribu- 
tion P({A„}). Within very general assumptions, such 
distribution was shown to satisfy the DMPK equation"^^ 
for disordered wiresi^ In the next section, we show that 
this also applies to graphene with ID disorder, and we 
use the known solution to such ID DMPK equation to 
compute linear response transport properties. 

Given the probabihty distribution P({A„}) for the 
transmission probabilities, a particularly useful object is 
the average density distribution p{X) = ^ni^i^ ~ ^n)) 
of the A„ above. It allows to compute any current fluc- 
tuation average in linear response (C) = TrC(Ft) by 



(C) = / dXp{X)C (scch^X) 



(6) 



Two particularly interesting cases are the conductivity, 
for which C(T) = {ie^ /h){L/W)T, and the shot noise, 
C{T) = {ie^\V\/h)T{l - T). In this work, we will also 
consider the Fano factor, defined as the ratio of disorder 
averages F = (T{1 - T))/{T). 



III. ONE-DIMENSIONAL DISORDER 

We will first use the preceding transfer matrix for- 
malism to study the case of one-dimensional disorder 
V{x,y) = V{x). In this case, Eq. (|4]) does not mix 
modes, which will be labeled by a well defined transver- 
sal momentum q. The corresponding A, of Eq. ([5]) will 
depend on the realization of disorder but will be other- 
wise mutually independent. If V{x) is further modeled 
by piecewise constant potential steps of width Axi [see 
Fig. [1] (a)] , it is straightforward to compute the resulting 



We first use this approach to numerically compute the 
reduced probability distribution P{Xq) for a given mode 
q. Random potential realizations are sampled with fixed 
Axi — f (representing the disorder correlation length) 
and statistically independent Gaussian potentials V{xi) 
of variance Var[F(a;i)] = ay. Each of these independent 
potential steps models a ID version of a charge puddle 
in graphene. '^^ The histogram of the computed transmis- 
sions for each realization and mode gives the distribution 

The resulting distribution P{Xq) for samples of length 
i ^ C (large number of puddles) is rather simple, see 
inset of Fig. [2l For high transverse momenta q\/^L ^ 1, 
it evolves into a Gaussian centered at a large A, while 
at lower momenta, it evolves into a Rayleigh-type distri- 
bution. Such distribution makes the limit of large num- 
ber of puddles non-self-averaging, i.e., at low momenta 
(C(A,))^C((A,)). 

Analytically, one can actually prove that, in the limit 
of large samples L ^ ^, the distribution P{Xq) satisfies 
the single channel DMPK equation for graphene (see Ap- 
pendix El 



dP Id f dP 



for a given mode q. In general, this equation involves 
a certain (possibly q dependent) mean free path I to be 
determined. Its solution is known and is given by^SiiM 

r..,x „ /2"e-^/4 r°° ^ Me~"'/^sinh2A 



-y/ cosh 2u ~ cosh 2 A ' 



where s — s[q) ^ L/l(q). As a side note, this solution is, 
for most practical purposes, quite indistinguishable from 
the simpler Rice distribution, Ps{X) « y^I/i^'^)' 
fined as 



■ exp 



A2 + A2\ fXX 



(9) 



transfer matrix for a given q: Tq{L) = e 



Mqq(xi)Axi 



The DMPK equation does not describe transport 
through evanescent modes, which are explicitly neglected 
in its derivation, see Eq. (jA7p . Transport through clean 
graphene at zero doping, on the other hand, is domi- 
nated by evanescent modes. The weak disorder limit of 
the DMPK solution is therefore quite different from the 
case of clean graphene, which has P{Xq) ~ 5{Xq — qL)^ 
In fact, the validity of the DMPK solution requires long 
enough samples so that the high q evanescent modes 
that are unaffected by disor der have died out, while 
the smaller \q\ < -^/e^ + f modes are converted into 
propagating modes by the effect of disorder. The ig- 
nored evanescent modes have x-wave vector l/cj,! > cTyS,, 
so that, more quantitatively, this long sample condition 
reads UyLS^:^ 1. The weak disorder limit of the DMPK 
solution corresponds, therefore, to that in which a single 
charge puddle has a weak effect on the scattering matrix, 
(Ty^ ^ 1, but all of them together have a strong effect 
alLi > 1. 
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One way to find the correct form of s{q) in Eq. ^ 
is to directly derive the DMPK equation from Eq. 
This is done in Appendix |A] in the case of weak disorder 
and small doping, (Tv^,eS^ <C 1. A more powerful way, 
also applicable to strong disorder, cry^ ^1, is to obtain 
an exact result for some expectation values of a function 
of Ag, such as (cosh2Ag). It is possible to compute this 
average analytically from Eq. ^ in the case of ID piece- 
wise disorder with a large number of charge puddles, see 
Appendix |B] By comparing it to the same average ob- 
tained from the distribution [Eq. ([S])], one arrives at the 
relation 



for L ^ ^ and ^ 1, where 



1 



(10) 



(11) 



if CTy^ ^ 1, (Jy^L ^ 1 (we call this 'weak disorder'), and 
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{avO''i2V2TTavC + tt - 2) - {e^YiiT + V2navQ 



if "'vC ^ 2 (strong disorder). All the dependence on dop- 
ing and disorder details is contained in 7. The derivation 
of the above expressions for 7 is valid for piecewise gaus- 
sian disorder to second order in the average doping e and 
L ^ Note that an identical result is non-trivially ob- 
tained if a different average such as (cosh^ 2Aq) is used for 
the calculation, which indicates that Eq. ([8|) is not merely 
an approximation, but is indeed exact for small dopings 
<C 1, both in the strong and weak disorder limits. We 
conjecture that any other model of ID smooth-disorder 
would satisfy, for large enough samples and low dopings, 
the above power-law equation for s(q) cx with a cer- 
tain 7 (independent of q and L) , regardless of the details 
of disorder. 

Equation ([5]) is valid for any ratio W/L. In the case 
of wider-than-long graphcne sheets, however, the average 
eigenvalue density function Pid(A) = '^qPqW takes on 
a remarkably simple form 



Pid(A) 
Pid(A) 




Pid(A), 



V2 



Pid(0) 



Ki{u) sinh2A 
V cosh 2u — cosh 2 A 



(13) 



du, 



valid in practice for W > L. Here, Ki{u) is the first mod- 
ified Bessel function of the second kind. Note that for 
most practical purposes, the simpler function Pid(A) w 
e~^/o(A), where Iq is the modified Bessel function of the 
first kind, can be used to within excellent accuracy. Im- 
plicit in the derivation of Eq. (|13p is that all modes 
— 00 < g < 00 are properly described by Eq. which 
as discussed above is valid if <Jy^L ^ 1. 
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FIG. 3: (Color online) Numerical results for the average con- 
ductivity vs. system length L/^ through a graphene sheet 
with ID disorder. (Note that data with L > W are not 
shown.) The results are normalized to the ballistic minimal 
conductivity o"™'" = Ae^ /nh. A a x \fL]\ scaling is obtained 
already from L « ^. At low dopings, Eq. p4p -black lines- 
agrees with the numerical results. 
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FIG. 4: (Color online) Fano factor vs sheet length L for 
various system parameters and ID disorder. The saturation 
value is universal at low dopings (in red) and lies at the value 
F = 0.243 derived from Eq. (fT3)) . 



In Fig. [21 we compare the numerical Pid(A) as ob- 
tained by disorder sampling in ID to the above analyti- 
cal expression [Eq. I|13p ] without any fitting parameters, 
finding a very good agreement, even when W is not much 
bigger than L. 

One of the main features of piD for graphene with 
smooth disorder [Eq. (fT5|) ] is the absence of a localization 



length. In other words, conductivity defined as a 
grows indefinitely with system size, scaling as 



G 



w 




(14) 



for large i/^. A similar scaling in one dimension was 
found in the case of white noise disorder within an ap- 
proximate self-averaging assumptioufSS This is in con- 
trast to the constant a of ballistic graphene-^- and to the 
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a ^ logL scaling behavior observed for 2D disorder by 
Bardarson et al.,^^ as we will discuss in the next section. 
The L-^^^ scaling of conductivity is not surprising, how- 
ever, in view of the property l{q) cx q~^- Indeed, within 
the window of propagating q modes, those with sizable 
transmission (ballistic modes) satisfy l{q) ~ l/(,q^ }^ L. 
Therefore, there will be approximately W/ ^/^L such bal- 
listic modes (those with smaller q). These modes will 
dominate transport, so that one expects the same L 
scaling for the conductance and hence a conductivity 
a oc \/L. 

Another important feature of Eq. (fTS]) is that all the 
details of the system (disorder strength, size, and dop- 
ing) enter as a A-independent prefactor to the density 
/5id(A) and do not affect the shape of the density pro- 
file. This directly implies that, in the parameter regime 
of validity for Eq. (|13p . the ratio of any pair of cur- 
rent cumulants is universal (independent of system size, 
mean free path, and disorder strength) since all current 
cumulants will scale with these parameters in the same 
way as the conductivity. In particular, the Fano factor 
in the presence of ID puddle disorder close to the Dirac 
point becomes Fid ~ 0.243, below the pseudodiffusive 
prediction F — 1/3 for ballistic graphenej^ The numeri- 
cal results for the conductivity and Fano factor in ID are 
shown in Figs. [3]and|31 respectively, as a function of sys- 
tem length L/S^. The numerical Fano factor is indeed seen 
to saturate close to F = 0.243 independent of system pa- 
rameters. The average conductivity also scales as \/~L]\ 
at low dopings as expected. Interestingly, this scaling 
persists also at higher dopings for which the derivation 
of s((j) [Eq. pU]) ] ceases to be valid. 



IV. TWO-DIMENSIONAL DISORDER 

In a realistic model of disordered graphene with W > 
L, mode mixing becomes important, so the preceding 
discussion of ID disorder need not apply. Indeed, as 
mentioned in the Introduction, some authors claimed^^ 
that large graphene sheets would, in fact, reach a uni- 
versal value of conductivity in the presence of smooth 
disorder (no intervalley scattering), while others^ nu- 
merically found a nonuniversal conductivity scaling as 
a oc log L for deltalike and Gaussian-correlated 2D disor- 
der potentials. As we will show in this section, our results 
support the latter nonuniversal conductivity in puddle 
disordered graphene. In addition, our results suggest the 
existence of universal ratios of current cumulants in the 
limit W > L ^ S^, just as in the ID case. 

In the presence of mode mixing due to two-dimensional 
disorder, the transfer matrix technique can still be used 
for numerical simulations, although it becomes more 
computationally intensive. Our goal once more is to com- 
pute the density P2d(A) = J2n^nW presence of 
charge puddle disorder of typical size ^ in both x and y 
directions, which we model by a random potential profile 
V{x,y) = '^gCiq{x)e^''^ . Harmonics aq{x) are approxi- 
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FIG. 5: (Color online) Normalized transmission density dis- 
tribution for 2D disorder and zero doping. For large enough 
systems, all the curves coincide (within the noise), suggest- 
ing a universal scaling for current cumulants just as for ID 
disorder (red dotted line). 

mated by piecewise constant functions in x at intervals 
of size ^. For a given x, aq{x) is assumed to be Gaus- 
sian distributed with variance Var[ag(a;)] = ay for all 
q < 27r/^. Higher harmonics are suppressed, as expected 
by the smoothness of charge puddles of typical size ^. The 
resulting potential, depicted in Fig. [Hb), is convenient 
for numerical computations [Vqiq{x) — aq'-q{x))] but can 
still be considered realistic for charge puddle disorder. 
For practical calculations, a high-momentum cutoff must 
be introduced; it is chosen high enough so that the re- 
sult for P2d(A) is cutoff independent. The computation 
of the transmissions is performed by composition of scat- 
tering matrices, rather than by multiplication of transfer 
matrices, since the latter method is unstable when many 
modes (higher momenta) have small transmissions, as is 
the case here. 

As can be seen in Fig. [51 the density /02d(A) clearly 
deviates from the ID solution (red dotted line). Strik- 
ingly, however, the shape of /32d(A), up to a global rescal- 
ing factor A2D (see below), is clearly independent of the 
sample length, width, or scattering strength as long as 
L/^ ^ 1 (complete mode mixing) and W > L (con- 
tinuum of transport momenta). When such conditions 
are satisfied, the disordered system acquires some uni- 
versal features. All the dependence of /92d(A) on system 
parameters W^L,^, ay, and e appears to enter as a A- 
independent prefactor, 

P2d(A) « A2D{L,W,^,av,e)p2DiX), (15) 

for some pure function of A, P2bWi just as in the ID case 
Eq. (jl3p . It is interesting to compare this pure function, 
plotted in Fig. O to the Dorokhov result'^^''^^ for diffusive 
metals which is A independent for large samples. Ballistic 
graphene (zero disorder) also has a A-independent p(A). 

Conductance and higher current cumulants are ob- 
tained from Eq. In contrast to the ID case, we 
see in Fig. [H] that the minimal conductivity GL/W is 
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FIG. 6: (Color online) Conductivity in 2D disordered 
graphene for various average dopings e and system param- 
eters. A clear logL scaling is observed up to L ~ where 
the transition to a quasi- ID ribbon regime begins. 
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FIG. 7: (Color online) As for ID disorder, the Fano factor 
saturates to a system independent value, this time around 
F = 0.295 for large systems sizes. 



proportional to logL/^, in agreement with Ref. [2J. The 
observed logarithmic dependence of the conductance on 
the sample size L might be attributed to the weak an- 
tilocalization. According to Ref. (see also Ref. |28| ) , the 
weak antilocalization correction to the conductance of a 
2D disordered graphene in the diffusive regime L 3> 
with no intervalley scattering, no magnetic impurities, 
and no interaction between electrons, should read 



da: 



2e2 , L 
— In - = CT 



r"in 



mean free path 1(e) should grow with energy as it did in 
case of ID disorder. 

From the logarithmic dependence of the conductiv- 
ity, it immediately follows that A2r,{L,W,£,,av,£) oc 
iog(£/$) ^ rpjjjg j]-Qpiies ilig^i shot noise and higher current 
cumulants should all scale with '"s^^/'^) . This behavior is 
indeed observed numerically (not shown explicitly), con- 
firming Eq. (|15p . 

Just as in the ID case, the form of Eq. (fT5|) implies 
once more that the ratio of any pair of current cumulants 
will be universal if L/^ 1 and W > L. The Fano 
factor, in particular, acquires a universal value i^2D ~ 
0.295, as can be seen in Fig. [T] Note that 1/3 > F2D > 



V. CONCLUSIONS 

In conclusion, we have characterized the full transport 
statistics of disorder graphene in the absence of inter- 
valley scattering. We have computed numerically the 
transmission probability distribution both for ID and 2D 
disorders. In both the ID and 2D cases we have found 
non-self-averaging statistics and no localization length. 
All current cumulants are seen to grow monotonously 
with system length without any saturation, which sug- 
gests that experimental results for the minimal conduc- 
tivity should be limited by contact resistance. 

In the ID case, we have shown that the solution for 
long samples satisfies a single channel DMPK equation, 
for which we give the exact analytical solution and the 
expression for the mean free path at low dopings. Such 
mean free path scales as I oc with the transverse mo- 
mentum, which has important consequences for the scal- 
ing properties of current cumulants. Such scaling prop- 
erties, observed numerically also in 2D, suggest that, for 
long and wide enough samples, all cumulants scale in an 
identical fashion with system size, leading to universal 
ratios of any two cumulants. Although the same phe- 
nomenon is obtained in the theory of diffusive metals;^ 
the resulting transport statistics is less trivial in the case 
of graphene. The Fano factor in 2D disordered graphene 
F = 0.295, for example, is below the value in diffusive 
metals and clean graphene F — 1/3. 



(16) 



Our low energy numerical data appear to be in reason- 
able agreement with this expression. At < 4, we find 
da/dlnL ^ 1.3 - 1.5 a^'" (see Fig. 6). Unfortunately, 
we do not know the value of the mean free path I and 
therefore cannot independently verify if our samples are 
in the diffusive regime. Most probably, at low energies, 
we get I < L and the diffusive model applies, while at 
higher energies, we have I ^ L which leads to a stronger 
dependence of the conductance on L. Such an interpre- 
tation of the results plotted in Fig. [S] suggests that the 
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APPENDIX A: DERIVATION OF THE 
ONE-DIMENSIONAL 
DOROKHOV-MELLO-PEREYRA-KUMAR 
EQUATION FOR GRAPHENE 

Following Titov, — we parametrize the T— matrix as 
follows: 



r = 



e^'^icoshAo e"*'^^sinhA 



From Eq. we derive the following equations^^ 

^ = qcos26q, 
d9 

= e-V(x) - qsm26„ coth2A„, 

ax 



(Al) 

(A2) 
(A3) 



where 9q — {ipi — ip2)/2. The initial condition for Eqs. 
(E2)) and jM]) reads \q{0) = and 6*5(0) = 0. Our aim is 
to derive an equation for the averaged distribution func- 
tion P{Xq, 9q). We proceed along the standard route and 
begin with the equation for the non-averaged distribu- 
tion function ^(Ag, 0,) — 5[\q — \q{x)]5[9q — 0q{x)\, whcrC 

Ag(x), 9q{x) are the solution of Eqs. and This 
equation reads 



dP 

dx 



d 

dfq 
d 



((e - V{x) - q sin 261, coth 2Xq)P) 



(q cos 29q P) . 



(A4) 



Now we assume that q£,,e£,,crv£, ^ 1- In this case, the 
fluctuating potential can be treated as d correlated, i.e., 
we consider 



{Vix,)V{x2)) ^ a^^S{xi - X2). 



The averaging over V{x) in Eq. (|A4[) becomes very sim- 
ple and we arrive at the usual Fokker-Plank equation for 
P=(P): 



dP 
dx 



dP crUd^P 



'd9q 2 



891 



q cos 29q 



dP 



+ qcoth.2\q 
Next, we assume 



2 2 



2cos29qP- 



sin26', 



dP 



^4 c2 



(A6) 



(A7) 



Then, the isotropization of P{Xq,9q) over the angle 9q 
happens fast, so that derivatives with respect to 9q in 
Eq. (jA6p are small. We split the distribution P{Xq,9q) 
into the sum of an isotropic part and small anisotropic 
correction, 

Pi\q,9q) = P{\q) + a{Xq) COS 29q + /3(A,) siu 29q. (A8) 

We derive three coupled equations for P, a and (3. To 
this end, we first average Eq. (jA6p over the angle 9q, 



then multiply it with cos 29q and average the result, and, 

with sin 29q and average 



finally, we multiply Eq. I|A6 
over 9q again. We then get 



dP 

dx 
1 da 
2'dx 

2'dx 



q da 
2d)Z 



q dP 



-el3 - ay^a - ^ ^ + q coth 2A, P, 



(A9) 



Under the condition [Eq. (|A7p ] we can set da/dx — 
df3/dx = 0. Then a and /3 are easily excluded and we 
arrive at the one channel DMPK equation 



dP Id f dP 
with the effective mean free path 



I 



(AlO) 



(All) 



Let us briefly discuss this result. First of all Eq. (All) 
predicts an infinite mean free path for an electron mov- 
ing perpendicular to the potential barriers [q — 0), i.e. 
no back-scattering occurs in this case. This is a man- 
ifestation of Klein paradox. At finite values of q the 
momentum of an incident electron is no longer perpen- 
dicular to the surface of the barrier, which makes the 
back-scattering possible. As a result the mean free path 
(All) becomes finite. Such behavior of the mean free 
path is encoded in the mathematical structure of Eqs. 
([2]) and ([3]). Namely, it is related to the fact that, in 
the case of ID disorder, the off-diagonal matrix elements 
(A5) oi Mqiq, responsible for back-scattering are proportional 



'99 

to q. One can actually derive the mean free path (All) 
directly from Eqs. ©-(HI in a simple way. To this end 
one should assume q to be sufficiently small and treat the 
off-diagonal elements of the matrix Mqiq perturbatively. 
Since Eq. ([2]) is formally similar to the Schrodinger equa- 
tion for a spin rotating in magnetic field, the "rate" of 
back-scattering, which should be identified with 1/21 in 
our problem, is given by a Fermi golden rule like expres- 
sion 



1 
2/ 



i^Re 



dy 



2^e(x-y)-2^f^ dzV{z) 



V{x) 



.(A12) 



The averaging over the fluctuating potential V{x) with 
the correlator (jASp reduces to the evaluation of a simple 
gaussian path integral and indeed leads to the expression 
(jAll|l for the mean free path. 

Note that the preceding derivation requires white- 
noise-like disorder with ayS, <C 1. We state without proof 
that, in the limit of strong disorder ay£, > 2, Eq. (|A10[) 
still holds, as is made plausible by the consistent deriva- 
tion of s = L// in Appendix iBl In this case however, the 
mean free path has a different form from Eq. (jAlip . 
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Moreover, the condition (jA7p impHes that the DMPK 
equation apphes only to modes which lie within a window 



< 



All modes with higher q are weakly 



sensitive to disorder and are, in fact, evanescent. When 
this window is finite (for finite disorder or finite doping) , 
the evanescent modes can be ignored provided that the 
sample is long enough. The precise condition is OyL^ ^ 
1. In this limit, the contribution of the evanescent modes 
with |gr| > a/c^ + CTy^^ to the conductance is negligible 
and Eq. becomes valid. 



APPENDIX B: ANALYTICAL RESULTS FOR 
ONE-DIMENSIONAL DISORDER AVERAGING 

The problem of computing the mean free path l{q) that 
enters the ID distribution [Eq. ^] is tackled here by 
finding analytical solutions for the expectation value of 
a certain function of A, in our case (cosh 2 A), and then 
adjusting s = L/l in order to recover that same result 
from the distribution ([5]). This is done for low dopings 
and both for the strong and the weak disorder limits. 
The procedure is carried out also for a different func- 
tion (cosh^ 2A), which yields an identical mean free path, 
which confirms the fact that Eq. ([5]) is the exact distribu- 
tion for transmissions through graphene with ID smooth 
disorder, as modeled b y E q. 

As explained in Ref. [26|, an alternative form of Eq. ([2]) 
for mode q under ID disorder is given by Eqs. (|A2p and 
(IA3p . The dynamics of 9 and A therein is coupled. How- 
ever, we are interested in expectation values of functions 
that do not involve 6 (as is the case of any observable 
current cumulant). The following exact reformulation of 
Eqs. (jA2p and (|A3[) proves useful to obtain them: 



v{x) = 

Q{x) = 

dv{x) 
dx 



cosh2A 
sinh 2A cos 29 
sinh 2A sin 29 

2q 
2q 

2{e-V{x)) 

Q{x)v{x), 



(Bl) 







-2{e-V{x)) I ,(B2) 


(B3) 



with the exact solution 



v{L) = Pcxp 



Q{x)dx v(Q). 



(B4) 



where P exp stands for a path ordered exponential. With 
a piecewise potential V{x) in which each step of size ^ is 
statistically independent of the others, the average {v{L)) 
reads 



(«(L)) = [(eQ«)]'^/«t;(0) 



(B5) 



We therefore need to obtain the average oiU = e^^ for a 
single charge puddle. To proceed analytically we have to 



assume that the average doping e is small, so we expand 
U to second order in efSi We will furthermore assume 
that the system contains a large number of puddles, N = 
L/^ > 1. In this will become apparent later, it 

is also sufficient to expand U to second order in q, since 
higher orders will not contribute to (U)^ in the large N 
limit. 

Let us now consider the case of strong disorder, which 
as it turns out means ayS. ^ 2. After the two previous 
expansions, U can be averaged by using a Gaussian dis- 
tribution for V with dispersion cry . The resulting expres- 
sion contains terms of the form Erf (v^cry^) and e"^"'^^ , 
which can be greatly simplified to 1 and respectively 
in the strong disorder limit ayS. ^ 2. To exponentiate 
the result one computes the eigenvalues of the resulting 
([/). One eigenvalue is zero, another is small as q^ and 
the last one is close to 1, so that if D is the matrix that 
diagonalizes (C/), one can write 



[(e««)] 



(1 + (gO'i)"^/^ 

= D\ M?!]"-^^ 





where 
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D- 



(B6) 



{avO^{2V2^av^ + tt - 2) - (eg)^(7r + V2^avO 



n{eO^-{avOHn^2) 
2(ayO' 



In the limit of large number of puddles N — L/£^ ^ oo, 
we have [l + iqO^l]^^^ ^ exp[q^^Lj], and [{qO^l]^^^ 
0, and higher order corrections in q become irrelevant. 
By writing down the expression for D in leading (zeroth) 
order in q, and selecting the first element of {v{L)) in Eq. 
(IB5I). we arrive at 



(cosh 2 A) = exp [q^L^^) 



(B7) 



From this result it becomes clear that is indeed a small 
expansion parameter for the relevant momenta, since the 
typical momentum scale that appears after averaging is 

<z^i/vxe«i/e 

Similarly we can obtain (cosh 2A) by doing an analo- 
gous computation for V2{L) = v{L) ® v{L), whose equa- 
tion of motion involves Q (g) 1 -I- 1 (g) Q instead of Q. The 
result in that case for the first element of {v2{L)) reads 



1 



(cosh^ 2A) = - + - exp (3g2L^7) 



(B8) 



The above is valid for strong disorder. In the case of 
weak disorder, we can expand 



exp 



|iog(^i + (Q)e + i(g^)e' 



exp<; I 
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Once again, only one of the eigenvalues of the exponent 
is relevant for the computation of (cosh 2A) . In the limit 
of large L/S^, this eigenvalue reads q'^^L'j, where now 



and we again find 

(cosh2A) 
(cosh^ 2A) 



exp {q'^L£,'-f) , 
1 2 



exp (Sq^L^^) 



(B9) 



The same averages computed from the ID DMPK dis- 
tribution [Eq. ([5])] read 



(cosh2A) = exp(2s) 
1 2 

(cosh^ 2A) = - + - exp(6s). 
3 3 



(BIO) 



Comparing Eqs. (|B7| - (|B10p . we get s = ^q'^L^j, i.e. Eq. 
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